library(lme4)
library(lmerTest)
library(afex)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(gt)
library(purrr)
library(stringr)
library(tidyr)
library(hms)
library(lubridate)
library(ggeffects)
mid_ps <- read.csv("../../../data/midway_2024/transformed/midway_ek_alpha_normalized_2024.csv")
mid_ps <- mid_ps %>%
rename(ID = specimenID)
mid_ps$date <- ymd(mid_ps$Date)
mid_ps$rlc_end_time <- as_hms(mid_ps$rlc_end_time)
mid_ps$run <- as.factor(mid_ps$run)
mid_ps$salinity <- as.factor(mid_ps$salinity)
mid_ps$nitrate <- as.factor(as.character(mid_ps$treatment))
mid_ps$deltaNPQ <- as.numeric(mid_ps$deltaNPQ)
mid_ps <- mid_ps %>%
group_by(day_group) %>%
mutate(rank = rank(rlc_end_time)) %>%
ungroup()
mid_ps <- mid_ps %>%
group_by(day_group) %>%
mutate(rlc_order = floor((rank -1)/3) +1)
mid_ps <- mid_ps %>%
mutate(treatment = case_when(
nitrate == 1 & salinity == 35 ~ "1",
nitrate == 1 & salinity == 28 ~ "2",
nitrate == 2 & salinity == 35 ~ "3",
nitrate == 2 & salinity == 28 ~ "4",
nitrate == 3 & salinity == 35 ~ "5",
nitrate == 3 & salinity == 28 ~ "6",
nitrate == 4 & salinity == 35 ~ "7",
nitrate == 4 & salinity == 28 ~ "8",
))
mid_ps_d9 <- subset(mid_ps, rlc_day == 9) #do not separate by plant_part
canopy <- subset(mid_ps, rlc_day == 9 & plant_part == "canopy")
canopy$treatment_graph[canopy$treatment == 1] <- "1) 0.5μmol/35 ppt"
canopy$treatment_graph[canopy$treatment == 2] <- "2) 0.5μmol/28 ppt"
canopy$treatment_graph[canopy$treatment == 3] <- "3) 2μmol/35 ppt"
canopy$treatment_graph[canopy$treatment == 4] <- "4) 2μmol/28 ppt"
canopy$treatment_graph[canopy$treatment == 5] <- "5) 4μmol/35 ppt"
canopy$treatment_graph[canopy$treatment == 6] <- "6) 4μmol/28 ppt"
canopy$treatment_graph[canopy$treatment == 7] <- "7) 8μmol/35 ppt"
canopy$treatment_graph[canopy$treatment == 8] <- "8) 8μmol/28 ppt"
under <- subset(mid_ps, rlc_day == 9 & plant_part == "under")
under$treatment_graph[under$treatment == 1] <- "1) 0.5μmol/35 ppt"
under$treatment_graph[under$treatment == 2] <- "2) 0.5μmol/28 ppt"
under$treatment_graph[under$treatment == 3] <- "3) 2μmol/35 ppt"
under$treatment_graph[under$treatment == 4] <- "4) 2μmol/28 ppt"
under$treatment_graph[under$treatment == 5] <- "5) 4μmol/35 ppt"
under$treatment_graph[under$treatment == 6] <- "6) 4μmol/28 ppt"
under$treatment_graph[under$treatment == 7] <- "7) 8μmol/35 ppt"
under$treatment_graph[under$treatment == 8] <- "8) 8μmol/28 ppt"
time_over_28 <- read_csv("../../../data/midway_2024/transformed/mins_28_plant_part.csv") #load dataset
mean_mins28 <- time_over_28 %>%
select(plant_part, run, mean_mins28_plant_part) %>% #keep only relevant columns of data
mutate(run = as.factor(run), mean_mins28_plant_part = as.factor(mean_mins28_plant_part))
canopy <- canopy %>%
left_join(mean_mins28, by = c("run", "plant_part")) #join the datasets
canopy <- canopy %>%
rename(mean_mins28 = mean_mins28_plant_part) #name is too long
under <- under %>%
left_join(mean_mins28, by = c("run", "plant_part")) #join the datasets
under <- under %>%
rename(mean_mins28 = mean_mins28_plant_part) #name is too long
growth_rate <- read.csv("../../../data/midway_2024/input/midway_growth.csv")
growth_rate$treatment <- as.factor(growth_rate$treatment)
growth_rate$d9_growth_percent <- (growth_rate$final_weight - growth_rate$initial_weight) / growth_rate$initial_weight * 100
d9 <- growth_rate %>%
select(ID, d9_growth_percent)#keep only relevant columns of data
canopy <- canopy %>%
left_join(d9, by = c("ID")) #join the datasets
under <- under %>%
left_join(d9, by = c("ID"))
check_model_fit <- function(model, terms) {
print(r.squaredGLMM(model))
hist(resid(model))
plot(resid(model) ~ fitted(model))
qqnorm(resid(model))
qqline(resid(model))
print(performance::check_model(model))
#print(tab_model(model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE,
#show.df = TRUE, show.zeroinf = TRUE))
plot(ggpredict(model, terms = terms))
}
get_data_means <- function(data, predictors_list, response) {
data %>%
group_by(across(all_of(c("nitrate", setdiff(predictors_list, "nitrate"))))) %>%
summarise(across(all_of(response), list(mean = mean), .group = "drop"))
}
compare_lmer_models <- function(data, response, predictors_list, random_effects) {
response <- as.character(response) # Ensure response is a string
models <- list() # Create a list to store the models
anova_result <- list() #initialize anova
random_effects_part <- paste0("(1 | ", random_effects, ")", collapse = " + ") #random effects used for all
fixed_effects_all <- paste(predictors_list, collapse = " + ") #make a formula for all predictors
formula <- as.formula(paste(response, "~", fixed_effects_all, "+", random_effects_part))
model_all <- lmer(formula = formula, data = data, REML = FALSE) #model_all is the version with both predictors for comparison
# Loop through predictors_list to create and fit models
for (i in seq_along(predictors_list)) {
fixed_effects <- paste(predictors_list[[i]], collapse = " + ")
formula <- as.formula(paste(response, "~", fixed_effects, "+", random_effects_part))
models[[i]] <- lmer(formula = formula, data = data, REML = FALSE) # Fit the lmer model and store in the list
anova_result[[i]] <- anova(models[[i]], model_all)
}
data_means <- get_data_means(data, predictors_list, response)
# Return results as a list
return(list(
models = models,
model_all = model_all,
anova_result = anova_result,
data_means = data_means
))
}
predictors_list <- c("salinity", "nitrate") #set the predictors list
canopy_pmax_results <- compare_lmer_models(
data = canopy,
response = "pmax",
predictors_list,
random_effects = c("plantID", "rlc_order")
)
summary(canopy_pmax_results$models[[1]]) # Model 1 summary (salinity)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 919.7 932.5 -454.8 909.7 91
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.79668 -0.55604 -0.09338 0.55720 2.60374
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 199.9 14.14
## rlc_order (Intercept) 157.1 12.53
## Residual 462.6 21.51
## Number of obs: 96, groups: plantID, 48; rlc_order, 32
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 101.765 4.821 50.311 21.107 <2e-16 ***
## salinity35 -7.549 6.116 39.693 -1.234 0.224
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## salinity35 -0.634
summary(canopy_pmax_results$models[[2]]) # Model 2 summary (nitrate)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 920.2 938.2 -453.1 906.2 89
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.71404 -0.57772 -0.06531 0.66096 2.58147
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 169.9 13.03
## rlc_order (Intercept) 147.3 12.14
## Residual 466.1 21.59
## Number of obs: 96, groups: plantID, 48; rlc_order, 32
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 109.257 6.530 50.666 16.732 <2e-16 ***
## nitrate2 -15.133 8.761 47.771 -1.727 0.0906 .
## nitrate3 -19.788 9.144 50.888 -2.164 0.0352 *
## nitrate4 -10.144 8.761 47.771 -1.158 0.2527
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) nitrt2 nitrt3
## nitrate2 -0.671
## nitrate3 -0.700 0.522
## nitrate4 -0.671 0.455 0.522
summary(canopy_pmax_results$model_all) # Model 3 summary (salinity and nitrate)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 920.6 941.1 -452.3 904.6 88
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.86204 -0.55814 -0.01586 0.65726 2.69009
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 151.4 12.31
## rlc_order (Intercept) 154.7 12.44
## Residual 463.5 21.53
## Number of obs: 96, groups: plantID, 48; rlc_order, 32
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 113.054 7.048 51.337 16.040 <2e-16 ***
## salinity35 -7.581 5.773 39.262 -1.313 0.1967
## nitrate2 -15.145 8.586 47.480 -1.764 0.0842 .
## nitrate3 -19.850 8.989 51.181 -2.208 0.0317 *
## nitrate4 -10.095 8.586 47.480 -1.176 0.2456
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.410
## nitrate2 -0.609 0.000
## nitrate3 -0.638 0.000 0.523
## nitrate4 -0.609 0.000 0.452 0.523
canopy_pmax_results$anova_result[[1]] # ANOVA did nitrate have an effect on Pmax?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 5 919.70 932.52 -454.85 909.70
## model_all 8 920.57 941.09 -452.29 904.57 5.1236 3 0.163
canopy_pmax_results$anova_result[[2]] # ANOVA did salinity have an effect on Pmax?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 7 920.25 938.20 -453.12 906.25
## model_all 8 920.57 941.09 -452.29 904.57 1.6766 1 0.1954
check_model_fit(canopy_pmax_results$model_all, terms = predictors_list)
## R2m R2c
## [1,] 0.08265697 0.44759
canopy_pmax_results$data_means
## # A tibble: 8 Ă— 3
## # Groups: nitrate [4]
## nitrate salinity pmax_mean
## <fct> <fct> <dbl>
## 1 1 28 108.
## 2 1 35 110.
## 3 2 28 101.
## 4 2 35 87.1
## 5 3 28 94.9
## 6 3 35 85.6
## 7 4 28 102.
## 8 4 35 95.0
under_pmax_results <- compare_lmer_models(
data = under,
response = "pmax",
predictors_list,
random_effects = c("mean_mins28", "rlc_order")
)
summary(under_pmax_results$models[[1]]) # Model 1 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 883.6 896.5 -436.8 873.6 91
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.39721 -0.58387 -0.02647 0.59074 2.70764
##
## Random effects:
## Groups Name Variance Std.Dev.
## rlc_order (Intercept) 15.89 3.986
## mean_mins28 (Intercept) 13.23 3.637
## Residual 501.46 22.393
## Number of obs: 96, groups: rlc_order, 16; mean_mins28, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 75.62709 4.26002 4.33413 17.753 3.26e-05 ***
## salinity35 0.05666 4.61214 90.54020 0.012 0.99
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## salinity35 -0.541
summary(under_pmax_results$models[[2]]) # Model 2 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 886.8 904.7 -436.4 872.8 89
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.37932 -0.63396 -0.05697 0.61647 2.84997
##
## Random effects:
## Groups Name Variance Std.Dev.
## rlc_order (Intercept) 12.06 3.473
## mean_mins28 (Intercept) 13.05 3.613
## Residual 500.01 22.361
## Number of obs: 96, groups: rlc_order, 16; mean_mins28, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 79.057 5.383 9.864 14.687 4.99e-08 ***
## nitrate2 -6.098 6.598 79.586 -0.924 0.358
## nitrate3 -4.421 6.694 45.972 -0.660 0.512
## nitrate4 -3.088 6.598 79.586 -0.468 0.641
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) nitrt2 nitrt3
## nitrate2 -0.613
## nitrate3 -0.622 0.507
## nitrate4 -0.613 0.485 0.507
summary(under_pmax_results$model_all) # Model 3 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 888.8 909.3 -436.4 872.8 88
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.37809 -0.63522 -0.05697 0.61709 2.85120
##
## Random effects:
## Groups Name Variance Std.Dev.
## rlc_order (Intercept) 12.06 3.472
## mean_mins28 (Intercept) 13.05 3.613
## Residual 500.01 22.361
## Number of obs: 96, groups: rlc_order, 16; mean_mins28, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 79.02837 5.85285 13.55302 13.503 3.03e-09 ***
## salinity35 0.05736 4.59678 90.32593 0.012 0.990
## nitrate2 -6.09755 6.59829 79.58591 -0.924 0.358
## nitrate3 -4.42113 6.69434 45.97127 -0.660 0.512
## nitrate4 -3.08787 6.59829 79.58591 -0.468 0.641
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.393
## nitrate2 -0.564 0.000
## nitrate3 -0.572 0.000 0.507
## nitrate4 -0.564 0.000 0.485 0.507
under_pmax_results$anova_result[[1]] # ANOVA did nitrate have an effect on Pmax?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 5 883.65 896.47 -436.82 873.65
## model_all 8 888.75 909.27 -436.38 872.75 0.8926 3 0.8272
under_pmax_results$anova_result[[2]] # ANOVA did salinity have an effect on Pmax?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 7 886.75 904.70 -436.38 872.75
## model_all 8 888.75 909.27 -436.38 872.75 2e-04 1 0.99
check_model_fit(under_pmax_results$model_all, terms = predictors_list)
## R2m R2c
## [1,] 0.009520906 0.05687984
under_pmax_results$data_means
## # A tibble: 8 Ă— 3
## # Groups: nitrate [4]
## nitrate salinity pmax_mean
## <fct> <fct> <dbl>
## 1 1 28 78.1
## 2 1 35 80.4
## 3 2 28 73.0
## 4 2 35 73.5
## 5 3 28 68.7
## 6 3 35 79.6
## 7 4 28 82.7
## 8 4 35 69.2
canopy_npqmax_results <- compare_lmer_models(
data = canopy,
response = "maxNPQ_Ypoint1",
predictors_list,
random_effects = c("plantID", "mean_mins28", "rlc_order")
)
summary(canopy_npqmax_results$models[[1]]) # Model 1 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 65.0 80.4 -26.5 53.0 90
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5570 -0.4536 -0.1271 0.4531 4.0118
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 0.042708 0.20666
## rlc_order (Intercept) 0.001894 0.04352
## mean_mins28 (Intercept) 0.012323 0.11101
## Residual 0.062045 0.24909
## Number of obs: 96, groups: plantID, 48; rlc_order, 32; mean_mins28, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.71940 0.09643 3.59648 7.460 0.0026 **
## salinity35 -0.03194 0.07855 46.01523 -0.407 0.6862
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## salinity35 -0.407
summary(canopy_npqmax_results$models[[2]]) # Model 2 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 42.4 62.9 -13.2 26.4 88
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1222 -0.4656 0.0082 0.4043 5.2251
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 0.0105110 0.10252
## rlc_order (Intercept) 0.0008604 0.02933
## mean_mins28 (Intercept) 0.0117902 0.10858
## Residual 0.0630749 0.25115
## Number of obs: 96, groups: plantID, 48; rlc_order, 32; mean_mins28, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.00567 0.09725 4.06386 10.341 0.000455 ***
## nitrate2 -0.36564 0.08413 46.87709 -4.346 7.40e-05 ***
## nitrate3 -0.37748 0.08442 41.65616 -4.472 5.87e-05 ***
## nitrate4 -0.46583 0.08413 46.87709 -5.537 1.35e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) nitrt2 nitrt3
## nitrate2 -0.433
## nitrate3 -0.434 0.502
## nitrate4 -0.433 0.497 0.502
summary(canopy_npqmax_results$model_all) # Model 3 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 44.1 67.1 -13.0 26.1 87
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1856 -0.4654 0.0101 0.3926 5.1965
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 0.0102212 0.1011
## rlc_order (Intercept) 0.0006758 0.0260
## mean_mins28 (Intercept) 0.0117925 0.1086
## Residual 0.0632712 0.2515
## Number of obs: 96, groups: plantID, 48; rlc_order, 32; mean_mins28, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.02164 0.10151 4.79409 10.064 0.00021 ***
## salinity35 -0.03210 0.05914 45.32554 -0.543 0.58987
## nitrate2 -0.36563 0.08386 46.79540 -4.360 7.08e-05 ***
## nitrate3 -0.37723 0.08408 41.28225 -4.487 5.68e-05 ***
## nitrate4 -0.46578 0.08386 46.79540 -5.555 1.28e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.291
## nitrate2 -0.413 0.000
## nitrate3 -0.414 0.000 0.501
## nitrate4 -0.413 0.000 0.497 0.501
canopy_npqmax_results$anova[[1]] # ANOVA did nitrate have an effect on NPQmax?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 6 64.984 80.371 -26.492 52.984
## model_all 9 44.068 67.147 -13.034 26.068 26.916 3 6.131e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
canopy_npqmax_results$anova[[2]] # ANOVA did salinity have an effect on NPQmax?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 8 42.362 62.876 -13.181 26.362
## model_all 9 44.068 67.147 -13.034 26.068 0.2933 1 0.5881
check_model_fit(canopy_npqmax_results$model_all, terms = predictors_list)
## R2m R2c
## [1,] 0.2745407 0.4660269
canopy_npqmax_results$data_means
## # A tibble: 8 Ă— 3
## # Groups: nitrate [4]
## nitrate salinity maxNPQ_Ypoint1_mean
## <fct> <fct> <dbl>
## 1 1 28 1.01
## 2 1 35 1.00
## 3 2 28 0.619
## 4 2 35 0.660
## 5 3 28 0.704
## 6 3 35 0.554
## 7 4 28 0.546
## 8 4 35 0.534
under_npqmax_results <- compare_lmer_models(
data = under,
response = "maxNPQ_Ypoint1",
predictors_list,
random_effects = c("plantID")
)
summary(under_npqmax_results$models[[1]]) # Model 1 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## -17.1 -6.8 12.5 -25.1 92
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.68725 -0.49727 0.06191 0.44741 2.92434
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 0.00591 0.07688
## Residual 0.03956 0.19889
## Number of obs: 96, groups: plantID, 48
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.57165 0.03272 48.00000 17.473 <2e-16 ***
## salinity35 0.08319 0.04627 48.00000 1.798 0.0785 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## salinity35 -0.707
summary(under_npqmax_results$models[[2]]) # Model 2 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## -35.3 -19.9 23.7 -47.3 90
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.09356 -0.48225 0.08922 0.55435 2.49521
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 0.00000 0.0000
## Residual 0.03577 0.1891
## Number of obs: 96, groups: plantID, 48
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.79508 0.03861 96.00000 20.595 < 2e-16 ***
## nitrate2 -0.25063 0.05460 96.00000 -4.590 1.34e-05 ***
## nitrate3 -0.21083 0.05460 96.00000 -3.862 0.000204 ***
## nitrate4 -0.26592 0.05460 96.00000 -4.871 4.38e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) nitrt2 nitrt3
## nitrate2 -0.707
## nitrate3 -0.707 0.500
## nitrate4 -0.707 0.500 0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
summary(under_npqmax_results$model_all) # Model 3 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## -38.1 -20.1 26.0 -52.1 89
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0301 -0.5251 0.1471 0.6169 2.7833
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 0.00000 0.0000
## Residual 0.03404 0.1845
## Number of obs: 96, groups: plantID, 48
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.75349 0.04211 96.00000 17.895 < 2e-16 ***
## salinity35 0.08319 0.03766 96.00000 2.209 0.029563 *
## nitrate2 -0.25063 0.05326 96.00000 -4.706 8.50e-06 ***
## nitrate3 -0.21083 0.05326 96.00000 -3.959 0.000145 ***
## nitrate4 -0.26592 0.05326 96.00000 -4.993 2.66e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.447
## nitrate2 -0.632 0.000
## nitrate3 -0.632 0.000 0.500
## nitrate4 -0.632 0.000 0.500 0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
under_npqmax_results$anova[[1]] # ANOVA did nitrate have an effect on NPQmax?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 4 -17.097 -6.8392 12.548 -25.097
## model_all 7 -38.065 -20.1147 26.033 -52.065 26.969 3 5.977e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
under_npqmax_results$anova[[2]] # ANOVA did salinity have an effect on NPQmax?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 6 -35.306 -19.920 23.653 -47.306
## model_all 7 -38.065 -20.115 26.033 -52.065 4.7592 1 0.02914 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_model_fit(under_npqmax_results$model_all, terms = predictors_list)
## R2m R2c
## [1,] 0.2808741 0.2808741
under_npqmax_results$data_means
## # A tibble: 8 Ă— 3
## # Groups: nitrate [4]
## nitrate salinity maxNPQ_Ypoint1_mean
## <fct> <fct> <dbl>
## 1 1 28 0.724
## 2 1 35 0.866
## 3 2 28 0.509
## 4 2 35 0.580
## 5 3 28 0.550
## 6 3 35 0.619
## 7 4 28 0.504
## 8 4 35 0.554
canopy_delta_npq_results <- compare_lmer_models(
data = canopy,
response = "deltaNPQ",
predictors_list,
random_effects = c("plantID")
)
summary(canopy_delta_npq_results$models[[1]]) # Model 1 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 47.9 58.1 -19.9 39.9 92
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6774 -0.5045 -0.1333 0.3854 4.1555
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 0.03378 0.1838
## Residual 0.06112 0.2472
## Number of obs: 96, groups: plantID, 48
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.53806 0.05178 48.00000 10.392 7.09e-14 ***
## salinity35 -0.02440 0.07322 48.00000 -0.333 0.74
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## salinity35 -0.707
summary(canopy_delta_npq_results$models[[2]]) # Model 2 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 24.4 39.8 -6.2 12.4 90
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6222 -0.5588 -0.1118 0.4499 5.3800
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 0.005754 0.07585
## Residual 0.061116 0.24722
## Number of obs: 96, groups: plantID, 48
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.81296 0.05501 48.00000 14.779 < 2e-16 ***
## nitrate2 -0.36958 0.07779 48.00000 -4.751 1.87e-05 ***
## nitrate3 -0.35367 0.07779 48.00000 -4.546 3.71e-05 ***
## nitrate4 -0.42513 0.07779 48.00000 -5.465 1.63e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) nitrt2 nitrt3
## nitrate2 -0.707
## nitrate3 -0.707 0.500
## nitrate4 -0.707 0.500 0.500
summary(canopy_delta_npq_results$model_all) # Model 3 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 26.2 44.2 -6.1 12.2 89
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5826 -0.5452 -0.0900 0.4251 5.3485
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 0.005605 0.07487
## Residual 0.061116 0.24722
## Number of obs: 96, groups: plantID, 48
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.82516 0.06138 48.00000 13.444 < 2e-16 ***
## salinity35 -0.02440 0.05490 48.00000 -0.444 0.659
## nitrate2 -0.36958 0.07763 48.00000 -4.761 1.81e-05 ***
## nitrate3 -0.35367 0.07763 48.00000 -4.556 3.60e-05 ***
## nitrate4 -0.42513 0.07763 48.00000 -5.476 1.57e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.447
## nitrate2 -0.632 0.000
## nitrate3 -0.632 0.000 0.500
## nitrate4 -0.632 0.000 0.500 0.500
canopy_delta_npq_results$anova[[1]] # ANOVA did nitrate have an effect on deltaNPQ?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 4 47.858 58.115 -19.929 39.858
## model_all 7 26.202 44.152 -6.101 12.202 27.656 3 4.29e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
canopy_delta_npq_results$anova[[2]] # ANOVA did salinity have an effect on deltaNPQ?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 6 24.399 39.785 -6.1995 12.399
## model_all 7 26.202 44.152 -6.1010 12.202 0.1971 1 0.6571
check_model_fit(canopy_delta_npq_results$model_all, terms = predictors_list)
## R2m R2c
## [1,] 0.3002212 0.3590055
canopy_delta_npq_results$data_means
## # A tibble: 8 Ă— 3
## # Groups: nitrate [4]
## nitrate salinity deltaNPQ_mean
## <fct> <fct> <dbl>
## 1 1 28 0.814
## 2 1 35 0.812
## 3 2 28 0.443
## 4 2 35 0.444
## 5 3 28 0.514
## 6 3 35 0.405
## 7 4 28 0.382
## 8 4 35 0.394
under_delta_npq_results <- compare_lmer_models(
data = under,
response = "deltaNPQ",
predictors_list,
random_effects = c("plantID")
)
summary(under_delta_npq_results$models[[1]]) # Model 1 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## -59.4 -49.1 33.7 -67.4 92
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.15735 -0.57743 -0.01169 0.42563 2.39882
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 0.01097 0.1047
## Residual 0.02006 0.1416
## Number of obs: 96, groups: plantID, 48
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.39404 0.02958 48.00000 13.32 <2e-16 ***
## salinity35 0.06402 0.04183 48.00000 1.53 0.132
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## salinity35 -0.707
summary(under_delta_npq_results$models[[2]]) # Model 2 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## -95.1 -79.7 53.5 -107.1 90
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1095 -0.4721 0.1215 0.4697 2.2747
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 0.0000 0.0000
## Residual 0.0192 0.1386
## Number of obs: 96, groups: plantID, 48
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.62083 0.02828 96.00000 21.952 < 2e-16 ***
## nitrate2 -0.26313 0.04000 96.00000 -6.579 2.49e-09 ***
## nitrate3 -0.23767 0.04000 96.00000 -5.942 4.51e-08 ***
## nitrate4 -0.27833 0.04000 96.00000 -6.959 4.21e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) nitrt2 nitrt3
## nitrate2 -0.707
## nitrate3 -0.707 0.500
## nitrate4 -0.707 0.500 0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
summary(under_delta_npq_results$model_all) # Model 3 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## -98.3 -80.4 56.2 -112.3 89
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.95853 -0.39173 0.03624 0.57564 2.10050
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 0.00000 0.0000
## Residual 0.01817 0.1348
## Number of obs: 96, groups: plantID, 48
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.58882 0.03076 96.00000 19.140 < 2e-16 ***
## salinity35 0.06402 0.02752 96.00000 2.327 0.0221 *
## nitrate2 -0.26313 0.03891 96.00000 -6.762 1.06e-09 ***
## nitrate3 -0.23767 0.03891 96.00000 -6.107 2.15e-08 ***
## nitrate4 -0.27833 0.03891 96.00000 -7.152 1.69e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.447
## nitrate2 -0.632 0.000
## nitrate3 -0.632 0.000 0.500
## nitrate4 -0.632 0.000 0.500 0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
under_delta_npq_results$anova[[1]] # ANOVA did nitrate have an effect on deltaNPQ?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 4 -59.359 -49.102 33.679 -67.359
## model_all 7 -98.319 -80.368 56.159 -112.319 44.96 3 9.437e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
under_delta_npq_results$anova[[2]] # ANOVA did salinity have an effect on deltaNPQ?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 6 -95.053 -79.667 53.526 -107.05
## model_all 7 -98.319 -80.368 56.159 -112.32 5.266 1 0.02175 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_model_fit(under_delta_npq_results$model_all, terms = predictors_list)
## R2m R2c
## [1,] 0.435658 0.435658
under_delta_npq_results$data_means
## # A tibble: 8 Ă— 3
## # Groups: nitrate [4]
## nitrate salinity deltaNPQ_mean
## <fct> <fct> <dbl>
## 1 1 28 0.543
## 2 1 35 0.699
## 3 2 28 0.351
## 4 2 35 0.365
## 5 3 28 0.350
## 6 3 35 0.417
## 7 4 28 0.333
## 8 4 35 0.352
canopy_ek_results <- compare_lmer_models(
data = canopy,
response = "ek.est",
predictors_list,
random_effects = c("plantID", "mean_mins28", "rlc_order")
)
summary(canopy_ek_results$models[[1]]) # Model 1 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 993.6 1009.0 -490.8 981.6 90
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8311 -0.5932 -0.0605 0.5266 3.3316
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 361.2 19.01
## rlc_order (Intercept) 240.9 15.52
## mean_mins28 (Intercept) 150.7 12.28
## Residual 1040.1 32.25
## Number of obs: 96, groups: plantID, 48; rlc_order, 32; mean_mins28, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 177.380 10.964 3.239 16.179 0.000329 ***
## salinity35 -13.284 8.712 39.775 -1.525 0.135214
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## salinity35 -0.397
summary(canopy_ek_results$models[[2]]) # Model 2 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 995.7 1016.2 -489.8 979.7 88
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1327 -0.5942 -0.0610 0.4520 3.0012
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 316.0 17.78
## rlc_order (Intercept) 258.3 16.07
## mean_mins28 (Intercept) 148.9 12.20
## Residual 1033.6 32.15
## Number of obs: 96, groups: plantID, 48; rlc_order, 32; mean_mins28, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 159.526 12.649 5.654 12.612 2.36e-05 ***
## nitrate2 12.348 12.498 46.880 0.988 0.3282
## nitrate3 7.310 12.981 50.099 0.563 0.5759
## nitrate4 25.190 12.498 46.880 2.016 0.0496 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) nitrt2 nitrt3
## nitrate2 -0.494
## nitrate3 -0.513 0.519
## nitrate4 -0.494 0.461 0.519
summary(canopy_ek_results$model_all) # Model 3 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 995.0 1018.1 -488.5 977.0 87
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.94436 -0.66754 0.00851 0.56066 3.12130
##
## Random effects:
## Groups Name Variance Std.Dev.
## plantID (Intercept) 247.8 15.74
## rlc_order (Intercept) 301.7 17.37
## mean_mins28 (Intercept) 145.2 12.05
## Residual 1020.1 31.94
## Number of obs: 96, groups: plantID, 48; rlc_order, 32; mean_mins28, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 166.014 13.063 6.446 12.709 8.39e-06 ***
## salinity35 -13.615 8.118 36.394 -1.677 0.1021
## nitrate2 12.618 12.067 45.994 1.046 0.3012
## nitrate3 7.773 12.626 50.791 0.616 0.5409
## nitrate4 25.734 12.067 45.994 2.133 0.0383 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.311
## nitrate2 -0.462 0.000
## nitrate3 -0.483 0.000 0.523
## nitrate4 -0.462 0.000 0.453 0.523
canopy_ek_results$anova[[1]] # ANOVA did nitrate have an effect on Ek?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 6 993.60 1009.0 -490.80 981.60
## model_all 9 995.04 1018.1 -488.52 977.04 4.5637 3 0.2067
canopy_ek_results$anova[[2]] # ANOVA did salinity have an effect on Ek?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 8 995.68 1016.2 -489.84 979.68
## model_all 9 995.04 1018.1 -488.52 977.04 2.6442 1 0.1039
check_model_fit(canopy_ek_results$model_all, terms = predictors_list)
## R2m R2c
## [1,] 0.07310292 0.448615
canopy_ek_results$data_means
## # A tibble: 8 Ă— 3
## # Groups: nitrate [4]
## nitrate salinity ek.est_mean
## <fct> <fct> <dbl>
## 1 1 28 158.
## 2 1 35 164.
## 3 2 28 178.
## 4 2 35 166.
## 5 3 28 176.
## 6 3 35 157.
## 7 4 28 196.
## 8 4 35 171.
under_ek_results <- compare_lmer_models(
data = under,
response = "ek.est",
predictors_list,
random_effects = c("mean_mins28", "rlc_order")
)
summary(under_ek_results$models[[1]]) # Model 1 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 970.8 983.7 -480.4 960.8 91
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9867 -0.6906 -0.1270 0.6353 3.0333
##
## Random effects:
## Groups Name Variance Std.Dev.
## rlc_order (Intercept) 36.832 6.069
## mean_mins28 (Intercept) 2.534 1.592
## Residual 1264.120 35.554
## Number of obs: 96, groups: rlc_order, 16; mean_mins28, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 133.936 5.489 6.557 24.402 1.1e-07 ***
## salinity35 -6.767 7.318 89.907 -0.925 0.358
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## salinity35 -0.667
summary(under_ek_results$models[[2]]) # Model 2 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 969.7 987.7 -477.9 955.7 89
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.12670 -0.68008 -0.08205 0.63288 2.99702
##
## Random effects:
## Groups Name Variance Std.Dev.
## rlc_order (Intercept) 99.844 9.992
## mean_mins28 (Intercept) 7.511 2.741
## Residual 1144.211 33.826
## Number of obs: 96, groups: rlc_order, 16; mean_mins28, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 114.119 7.986 17.418 14.291 4.68e-11 ***
## nitrate2 19.059 10.416 85.624 1.830 0.0708 .
## nitrate3 25.003 10.856 55.420 2.303 0.0250 *
## nitrate4 21.671 10.416 85.624 2.081 0.0405 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) nitrt2 nitrt3
## nitrate2 -0.652
## nitrate3 -0.680 0.521
## nitrate4 -0.652 0.457 0.521
summary(under_ek_results$model_all) # Model 3 summary
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: formula
## Data: data
##
## AIC BIC logLik deviance df.resid
## 970.8 991.3 -477.4 954.8 88
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.22521 -0.76232 -0.08919 0.57191 2.91970
##
## Random effects:
## Groups Name Variance Std.Dev.
## rlc_order (Intercept) 100.788 10.039
## mean_mins28 (Intercept) 7.803 2.793
## Residual 1131.695 33.641
## Number of obs: 96, groups: rlc_order, 16; mean_mins28, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 117.457 8.703 23.877 13.496 1.15e-12 ***
## salinity35 -6.748 7.004 86.873 -0.964 0.3379
## nitrate2 19.080 10.369 85.738 1.840 0.0692 .
## nitrate3 25.086 10.813 55.664 2.320 0.0240 *
## nitrate4 21.713 10.369 85.738 2.094 0.0392 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.402
## nitrate2 -0.596 0.000
## nitrate3 -0.621 0.000 0.521
## nitrate4 -0.596 0.000 0.456 0.521
under_ek_results$anova[[1]] # ANOVA did nitrate have an effect on Ek?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 5 970.83 983.65 -480.41 960.83
## model_all 8 970.79 991.30 -477.39 954.79 6.0408 3 0.1096
under_ek_results$anova[[2]] # ANOVA did salinity have an effect on Ek?
## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## models[[i]] 7 969.71 987.66 -477.86 955.71
## model_all 8 970.79 991.30 -477.39 954.79 0.9237 1 0.3365
check_model_fit(under_ek_results$model_all, terms = predictors_list)
## R2m R2c
## [1,] 0.07972884 0.1603012
under_ek_results$data_means
## # A tibble: 8 Ă— 3
## # Groups: nitrate [4]
## nitrate salinity ek.est_mean
## <fct> <fct> <dbl>
## 1 1 28 122.
## 2 1 35 111.
## 3 2 28 133.
## 4 2 35 136.
## 5 3 28 133.
## 6 3 35 139.
## 7 4 28 148.
## 8 4 35 123.
raw_plots <- function(data, response, label, pretty_color, aescolor, x, y, title, subtitle, ylim1, ylim2,
color1, color2, yint, vjust_t, hjust_t, vjust_s, hjust_s, labels1) {
histo <- data %>%
ggplot(aes(x = {{response}})) +
geom_histogram(fill = pretty_color, color = "black", size = 0.25, alpha = 0.85) +
labs(title = label) +
theme_bw()
plot <- data %>%
ggplot(aes(x = treatment_graph, y = {{response}})) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.75, size = 3, aes(color = {{aescolor}}), position = "jitter", show.legend = TRUE) +
labs(x = x, y = y, title = title, subtitle = subtitle) +
#scale_x_discrete(labels = c("0.5μmolN", "2μmolN", "4μmolN", "8μmolN")) +
ylim(ylim1, ylim2) + stat_mean() +
scale_color_manual(values = c(color1, color2)) +
scale_x_discrete(labels = labels1) +
geom_hline(yintercept=yint, color = "red", size = 0.5, alpha = 0.5) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = vjust_t, hjust = hjust_t),
plot.subtitle = element_text(face = "italic", size = 14, vjust = vjust_s, hjust = hjust_s))
return(list(
histo = histo,
plot = plot))
}
pmax_canopy <- raw_plots(
data = canopy,
response = pmax,
label = "Canopy",
pretty_color = "goldenrod1",
aescolor = salinity,
x = "Nitrate",
y = "Day 9 Pmax (μmols electrons m-2 s-1)",
title = "A - Canopy",
subtitle = "Chondria tumulosa",
ylim1 = 0,
ylim2 = 175,
color1 = "goldenrod3",
color2 = "darkgoldenrod1",
yint = 100,
vjust_t = -15,
hjust_t = 0.05,
vjust_s = -20,
hjust_s = 0.05,
labels1 = c("0.5 μmol", "0.5 μmol",
"2 μmol", "2 μmol",
"4 μmol", "4 μmol",
"8 μmol", "8 μmol")
)
plot(pmax_canopy$histo)
plot(pmax_canopy$plot)
pmax_under <- raw_plots(
data = under,
response = pmax,
label = "Understory",
pretty_color = "maroon4",
aescolor = salinity,
x = "Nitrate",
y = "Day 9 Pmax (μmols electrons m-2 s-1)",
title = "B - Understory",
subtitle = "Chondria tumulosa",
ylim1 = 0,
ylim2 = 175,
color1 = "deeppink4",
color2 = "deeppink3",
yint = 100,
vjust_t = -15,
hjust_t = 0.05,
vjust_s = -20,
hjust_s = 0.05,
labels1 = c("0.5 μmol", "0.5 μmol",
"2 μmol", "2 μmol",
"4 μmol", "4 μmol",
"8 μmol", "8 μmol")
)
plot(pmax_under$histo)
plot(pmax_under$plot)
npqmax_canopy <- raw_plots(
data = canopy,
response = maxNPQ_Ypoint1,
label = "Canopy",
pretty_color = "goldenrod2",
aescolor = salinity,
x = "Nitrate",
y = "Day 9 NPQmax",
title = "C - Canopy",
subtitle = "Chondria tumulosa",
ylim1 = 0,
ylim2 = 2,
color1 = "gold1",
color2 = "gold3",
yint = 1,
vjust_t = -15,
hjust_t = 0.05,
vjust_s = -20,
hjust_s = 0.05,
labels1 = c("0.5 μmol", "0.5 μmol",
"2 μmol", "2 μmol",
"4 μmol", "4 μmol",
"8 μmol", "8 μmol")
)
plot(npqmax_canopy$histo)
plot(npqmax_canopy$plot)
npqmax_under <- raw_plots(
data = under,
response = maxNPQ_Ypoint1,
label = "Understory",
pretty_color = "maroon1",
aescolor = salinity,
x = "Nitrate",
y = "Day 9 NPQmax",
title = "D - Understory",
subtitle = "Chondria tumulosa",
ylim1 = 0,
ylim2 = 2,
color1 = "hotpink2",
color2 = "hotpink4",
yint = 1,
vjust_t = -15,
hjust_t = 0.05,
vjust_s = -20,
hjust_s = 0.05,
labels1 = c("0.5 μmol", "0.5 μmol",
"2 μmol", "2 μmol",
"4 μmol", "4 μmol",
"8 μmol", "8 μmol")
)
plot(npqmax_under$histo)
plot(npqmax_under$plot)
delta_npq_canopy <- raw_plots(
data = canopy,
response = deltaNPQ,
label = "Canopy",
pretty_color = "goldenrod3",
aescolor = salinity,
x = "Nitrate",
y = "Day 9 deltaNPQ",
title = "E - Canopy",
subtitle = "Chondria tumulosa",
ylim1 = 0,
ylim2 = 2,
color1 = "lightgoldenrod2",
color2 = "lightgoldenrod3",
yint = 0.5,
vjust_t = -15,
hjust_t = 0.05,
vjust_s = -20,
hjust_s = 0.05,
labels1 = c("0.5 μmol", "0.5 μmol",
"2 μmol", "2 μmol",
"4 μmol", "4 μmol",
"8 μmol", "8 μmol")
)
plot(delta_npq_canopy$histo)
plot(delta_npq_canopy$plot)
delta_npq_under <- raw_plots(
data = under,
response = deltaNPQ,
label = "Understory",
pretty_color = "maroon3",
aescolor = salinity,
x = "Nitrate",
y = "Day 9 deltaNPQ",
title = "F - Understory",
subtitle = "Chondria tumulosa",
ylim1 = 0,
ylim2 = 2,
color1 = "hotpink1",
color2 = "hotpink3",
yint = 0.5,
vjust_t = -15,
hjust_t = 0.05,
vjust_s = -20,
hjust_s = 0.05,
labels1 = c("0.5 μmol", "0.5 μmol",
"2 μmol", "2 μmol",
"4 μmol", "4 μmol",
"8 μmol", "8 μmol")
)
plot(delta_npq_under$histo)
plot(delta_npq_under$plot)
ek_canopy <- raw_plots(
data = canopy,
response = ek.est,
label = "Canopy",
pretty_color = "goldenrod4",
aescolor = salinity,
x = "Nitrate",
y = "Day 9 Ek (μmols electrons m-2 s-1)",
title = "G - Canopy",
subtitle = "Chondria tumulosa",
ylim1 = 50,
ylim2 = 300,
color1 = "khaki4",
color2 = "khaki3",
yint = 150,
vjust_t = -15,
hjust_t = 0.05,
vjust_s = -20,
hjust_s = 0.05,
labels1 = c("0.5 μmol", "0.5 μmol",
"2 μmol", "2 μmol",
"4 μmol", "4 μmol",
"8 μmol", "8 μmol")
)
plot(ek_canopy$histo)
plot(ek_canopy$plot)
ek_under <- raw_plots(
data = under,
response = ek.est,
label = "Understory",
pretty_color = "maroon3",
aescolor = salinity,
x = "Nitrate",
y = "Day 9 Ek (μmols electrons m-2 s-1)",
title = "H - Understory",
subtitle = "Chondria tumulosa",
ylim1 = 50,
ylim2 = 300,
color1 = "palevioletred3",
color2 = "palevioletred4",
yint = 150,
vjust_t = -15,
hjust_t = 0.05,
vjust_s = -20,
hjust_s = 0.05,
labels1 = c("0.5 μmol", "0.5 μmol",
"2 μmol", "2 μmol",
"4 μmol", "4 μmol",
"8 μmol", "8 μmol")
)
plot(ek_under$histo)
plot(ek_under$plot)